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Abstract. The synchronization transition between two coupled replicas of spatio- 
temporal chaotic systems in 2+1 dimensions is studied as a phase transition into an 
absorbing state — the synchronized state. Confirming the scenario drawn in 1+1 
dimensional systems, the transition is found to belong to two different universality 
classes — Multiplicative Noise (MN) and Directed Percolation (DP) — depending on 
the linear or nonlinear character of damage spreading occurring in the coupled systems. 
By comparing coupled map lattice with two different stochastic models, accurate 
numerical estimates for MN in 2 + 1 dimensions are obtained. Finally, aiming to 
pave the way for future experimental studies, slightly non-identical replicas have been 
considered. It is shown that the presence of small differences between the dynamics of 
the two replicas acts as an external field in the context of absorbing phase transitions, 
and can be characterized in terms of a suitable critical exponent. 
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1. Introduction 

Synchronization of low dimensional chaotic oscillators, which have been intensively 
investigated in the last two decades [I], can be essentially regarded as a bifurcation 
from an uncorrelated to an entrained state. When single low dimensional oscillators 
are replaced by spatially extended systems exhibiting spatio-temporal chaos, however, 
synchronization becomes a genuine, fluctuation driven, phase transition which separates 
the uncorrelated and the completely synchronized states (2J |3]. The synchronization 
transition (ST) between two identical coupled replicas (starting from different, generic 
initial conditions) of the same spatio-temporal chaotic dynamics is thus a non- 
equilibrium critical phenomenon originating from deterministic dynamics. Chaotic 
fluctuations on the synchronized trajectory play the role of intrinsic stochastic terms, 
leading to diverging fluctuations as the critical point is approached [I]. 

Numerical studies in one dimensional dissipative extended systems have shown that 
the ST is essentially a phase transition into an absorbing state [5], i.e. the completely 
synchronized state, in which the two replicas of the extended system evolve on the same 
chaotic trajectory [21 E], similarly to low dimensional chaotic systems [6]. Noticeably, 
the observed continuous STs can be distinguished in two universality classes depending 
on the spatio-temporal propagation properties of the system close to criticality [2]. 
Consider a localized and finite perturbation (or synchronization error) to one of two 
otherwise identical (and thus synchronized) coupled replicas. At low coupling it will 
spread with a well defined propagation velocity vp implying a desynchronization 
of the two replicas. By definition, ST takes place at the critical coupling value for 
which vp vanishes, i.e. the synchronization error does not propagate anymore in space 
and time. On the other hand, the evolution of "infinitesimal" synchronization errors 
is well described by the linearized dynamics and can be characterized in terms of the 
asymptotic exponential growth rate, i.e. the so-called transverse Lyapunov exponent At 
[H,[6]: a positive (resp. negative) At implies exponential growth (resp. contraction) 
of infinitesimal perturbations. Whenever the local dynamics is sufficiently smooth, 
i.e. the linearization captures the essential features of the full dynamics, At vanishes 
together with vp at the critical coupling for the synchronization. When this happens ST 
belongs to the Multiplicative Noise universality class p] . When strong nonlinear 

effects make finite amplitude perturbations more unstable than innitesimal ones [10J, 
the transverse Lyapunov exponent vanishes before the error propagation velocity, so 
that at the ST transition one has At < and vp = 0. In this second case, which 
is closely related to the Stable Chaos phenomenon [TT] . ST belongs to the Directed 
Percolation (DP) universality class [5]. In the framework of one dimensional coupled 
map lattices (CMLs) - a common prototype of spatio-temporal chaotic systems - the 
critical behavior belongs to the DP universality class for (almost-) discontinuous local 
maps, and to the MN class for smoother maps [21 [3]. Remarkably, in the context of 



\ Or, to be more precise, to the so called MN1 class, see for instance Ref. jSJ. In this paper MN will 
always indicate the MN1 class. 
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CMLs models, it is possible to pass from DP to MN by varying a unique parameter, 
which controls the strength of nonlinear instabilities of the local map [T2] . 

Analytical arguments (coarse-graining techniques combined with linearization or 
finite-size analysis) show that these two STs can be mapped into the Langevin equations 
describing the MN or DP universality classes, respectively [3],[7j. It is worth noticing 
that these two universality classes have been originally identified in completely different 
contexts, such as epidemics spreading in reaction diffusion dynamics (DP) or the 
depinning of a fluctuating Kardar-Parisi-Zhang (KPZ) [13] interface (with a negative 
nonlinear term) from an underlying substrate (MN)1§|. Interestingly, while it is not 
possible to define a simple reaction diffusion system with a Markov dynamics belonging 
to the MN class, DP critical depinning can be observed when a short ranged attraction 
force is imposed between a KPZ interface and the underlying substrate [14J. Therefore, 
ST can be described in terms of a single Langevin equation [15]. Although it is yet 
unclear how the empirical Langevin equation for ST may eventually flow - as a control 
parameter is varied from the MN to the DP region - to the same DP fixed point as 
the Langevin equation for Reggeon field theory, these results strongly suggest that both 
universality classes may be described within a unique field-theoretic framework. 

In the last few years, ST has been object of intensive numerical investigations 
[21 El [71 [121 El HS1 EH EE] ; the analysis has also been extended to stochastic coupling 
[2], cellular automata [2)1 EI], systems with long-range interactions [22], [23] and delayed 
dynamical systems [2H |25], exploiting the analogy between the latter and spatially 
extended dynamics. So far, however, ST critical behavior has been investigated in 1 + 1 
dimensions only. In this respect, it should be also remarked that while DP critical 
exponents are known with great accuracy in any dimension, and MN exponents have 
been accurately measured in 1+1 dimensions [26] . the only known estimates of the 
MN critical indexes in 2+1 dimensions comes from direct numerical simulations of the 
Langevin dynamics [27]. Moreover, it is worth recalling that, for d < 2 the MN critical 
behavior is completely governed by a strong coupling fixed point [28] , so that no field 
theoretical estimates of critical exponents are currently available. 

The above situation is particularly unsatisfactory, especially considering that two 
dimensional spatio-temporal chaotic systems - such as chemical turbulence in quasi-two 
dimensional reactions [291 SO] or turbulent nematic liquid crystals |3lJ - are the most 
promising ones for studying ST in experiments. In this perspective, it is worth stressing 
that no experimental realization of the MN class has been so far realized, and it is only 
recently that the DP critical exponents have been measured in experiments: in 1+1 [32] 
and 2+1 dimensions [3T] . Therefore, besides the potential interest of experimentally 
realizing synchronization of chaotic extended systems, such an experiment would be the 
perfect ground for testing MN critical behavior in a physical framework. 

A major difficulty in devising a synchronization experiment is the assemblage of two 
perfectly identical replicas of the same system. In practice, it would be unavoidable to 

§ A second universality class, called MN2, can be denned starting from a KPZ with a positive nonlinear 
term. However, it is completely unrelated to synchronization problems. 
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experience very small differences, e.g., in the physical parameters entering the dynamics 
of the two systems. It is therefore of practical interest to quantify the influence of a small 
mismatch in the copies of the two systems. Another difficulty could lie in producing 
a suitable local coupling between the two replicas. However, this can be also realized 
through a stochastic forcing, for instance exposing two excitable chemical samples to 
the same random illumination (see, e.g., [30]), or by applying the same random external 
voltage to two replicas of the intermittent electrohydrodynamic convection regimes of 
Ref. |31j. As far as the critical properties of ST are concerned, deterministic or stochastic 
couplings are expected to share the same properties [2j [3] . 

This paper focuses on the synchronization transition in two spatial dimensions, 
within the framework of CMLs, with a twofold scope. First, we aim at verifying whether 
the scenario drawn in one dimensional systems for ST applies also to two dimensions. 
This requires accurate estimates of the MN critical exponents. For such reason we also 
investigate two stochastic models which are expected to belong to MN universality class. 
Second, mimicking what could happen in an experiment, we study ST in the presence 
of a small parameter mismatch between the two replicas. 

The paper is organized as follows. In Sect. El after presenting the CMLs models 
and recalling the definition and basic tools for measuring the critical exponents, we 
discuss the results of accurate numerical simulations for two different classes of maps, 
which in 1+1 dimensions have been proved to belong to the DP and MN universality 
classes. Sect. |3] focuses on accurate estimation of critical exponents for MN class in 
2 + 1 dimensions by means of two stochastic models. The section ends with a critical 
comparison between the computed exponents for the MN universality class and some 
scaling relations which have been put forward by previous studies. In SectJU mimicking 
what would typically be an experimental setting, we reconsider ST in the presence of 
a parameter mismatch and, accordingly, we introduce a new critical exponent for its 
characterization. Sect. [5] is devoted to final remarks and discussions. 

2. Synchronization transition in two spatial dimensions 

In the following we investigate the dynamics of two coupled replicas of two dimensional 
CMLs, defined on a L x L square lattice with periodic boundary conditions, evolving 
according to the dynamics: 

Ut+i(x, y) = (1 - -f)F(u t (x, y)) + jF(v t (x, y)) 

vt+i(x, y) = (1 - 7)F(v t (x, y)) + ^F(u t (x, y)) (1) 

with x, y — 1, 2, . . . , L. The variable 

zt(x, y) = V 2 £ z t (x, y) = (1 - e)z t (x, y) (2) 

+ -[z t (x +1,?/) + zt(x -l,y) + z t (x, y+l) + z t (x, y-l)], 

(with z — u, v) represents the nearest-neighbour diffusive coupling within each replica. 
Through this work, we set the diffusive coupling constant e equal to the democratic value 
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e = 4/5, that gives same weight to all neighbors. As from studies in the one dimensional 
version of Eq. ([1]), the DP or MN character of the ST relates to the functional form of 
the local map F(u) [21 E]- In particular, we have that if F(u) is continuous (such as the 
logistic or the tent map) the transition is in the MN universality class. Conversely, for 
discontinuous (or quasi-discontinuous [34]) maps, such as the generalized shift map 

F(u) = au mod 1 , (3) 

the ST belongs to the DP-class. The strength of the "transverse" coupling between the 
replicas is controlled by 7: for 7 = they are completely uncorrelated, while setting 
7 = 1/2 induces trivial complete synchronization in one time step, i.e. u t (x, y) = v t (x, y) 
for any t > 1 and for all x, y. Nontrivial ST, if present, is expected for a critical coupling 
value 7 C G]0 : l/2[ at which the synchronized state, u t (x,y) = v t (x,y), becomes stable 
(or, at least, marginally stable). For any 7 > j c , the two replicas (starting from different, 
generic initial conditions) converge toward the same spatio-temporal chaotic trajectory. 
In other terms, the synchronization error field w t (x,y) = \u t {x,y) — v t (x,y)\ tends 
towards zero for any x, y. The linear stability properties of the synchronized state 
are ruled by the transverse Lyapunov exponent At- For 7 > j c , At can be directly 
computed from the maximum Lyapunov exponent A of an uncoupled replica according 
to the relation [6] 

A T = A + ln(l - 2 7 ) . (4) 

The request At = thus determines the coupling 7 = (1 — exp(— A))/2, which coincides 
with the critical coupling for ST, i.e. 7 = 1c-, fo r maps belonging to the MN class. 
However, this result, based on the linear analysis, does not hold for systems belonging 
to the DP universality class. 

The suitable order parameter to characterize ST is the spatial average of the 
synchronization error w t (x,y) = \u t (x,y) — v t (x,y)\, i.e. 

1 L 

pit) = J2 Wt ( x >y) ■ ( 5 ) 

x,y=l 

Note that p(t) = if and only if w t (x,y) = for all x,y. If two replicas are identical 
(i.e. synchronized) at time to, they will remain so at all times t > to, implying that the 
synchronized state is absorbing, i.e. the dynamics cannot escape from it. 

Interestingly, the MN and DP absorbing states differ in their measure [33] which, 
for any finite system size L, is vanishing or finite, respectively. As a consequence any 
finite system belonging to the DP class falls into the absorbing state in a finite time. 
In CMLs, whose state variables are continuous, the synchronized state w t (x,y) = can 
only be reached asymptotically in time (barring computer round-offs), apparently at 
odd with the nature of the DP absorbing state. In Ref. [TJ, however, it has been shown 
that whenever the ST transition takes place at Vp = and At < (the DP case), 
all perturbations smaller than a certain finite (but vanishing for L — > 0) threshold are 
exponentially contracted towards zero, thus defining an effective finite measure for the 
absorbing state. 
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Figure 1. Snapshots of the synchronization error field w t (x, y) at four successive times 
(increasing from left to right) for a 2cHattice of 256 x 256 coupled Bernoulli shift maps. 
The two replicas are transversely coupled with 7 = 0.32824, slightly above the critical 
point. Color code goes from light toward black indicating decreasing intensity of the 
field (in logarithmic scale). 



We now consider the scaling theory for phase transitions into an absorbing state [5], 
which allows for defining the universal critical exponents, only depending on properties 
such as system symmetries and spatial dimensions, in terms of which we can characterize 
the ST. In the thermodynamic limit and at 7 = 7 C , one expects the order parameter to 
follow the asymptotic power law decay 

pit) ~ r e , (6) 

where 9 is a critical exponent. Close to the critical point, but slightly within the 
unsynchronized phase, i.e. for 7 < 7 C and 7 C — 7 <C 1, the asymptotic value of the 
order parameter p^ = linij^oo ^ J2t=i P$) ls characterized by the scaling relation 

Poo ~ (7c - if , (7) 

which defines a second critical exponent. Out of equilibrium critical phenomena are 
also characterized by the divergence of space and time correlations close to the critical 
point, in particular one expects 

£11 ~ It - TcP , £±~l7-7 c r x ( 8 ) 

where £11 (resp. is the temporal (resp. spatial) correlation length, where for simplicity 
we considered the system to be spatially isotropic. It can be shown [5] that only three 
exponents are independent, and in particular that the following relation holds 

9 = (3u\\ . (9) 

Furthermore, the ratio z = u»/v±, i.e. the so-called dynamical exponent, determines the 
relation between time and space correlations. The dynamical exponent can be measured 
exploiting the finite-size scaling relation 

p{t) = L- e *f{t/L z ) , (10) 

which holds at the critical point (/ being a universal scaling function). 
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Figure 2. Critical behaviors of ST for CMLs of Bernoulli maps, data refer to a 
1024 x 1024 lattice, (a) Order parameter as a function of time at criticality 7 C = 0.32817 
(black line) and slightly below (7 = 0.32815) and above (7 = 0.32819) the critical point 
(coloured thin lines). Data are averaged over 200 initial conditions. The dashed red 
line marks the t behavior with 9 — 0.449. (b) The compensated order parameter 
p(t)t e reveals a satisfactory plateau over about two decades, (c) Asymptotic order 
parameter p,^ as a function of the distance from the critical point. The dashed red 
line shows the best fit with (3 — 0.584(9). For each 7, is obtained averaging over 
100 initial conditions. 



2.1. Discontinuous Maps and Directed Percolation universality class 

We start analyzing the synchronization transition of two different replicas of the same 
CML (CQ) whose local dynamics is given by a discontinuous map. In particular, we 
consider the Bernoulli map ([3]) with a = 2. Fig. [1] displays the typical spatial 
structure of the difference field w t (x, y), as obtained iterating two replicas starting from 
(independent) random initial conditions, at successive times, for 7 slightly larger than 
the critical coupling value. The figure reveals typical DP patterns in proximity of the 
synchronization threshold. In order to make quantitative such observation, we need to 
measure the critical exponents. To this aim, we preliminarily determined the critical 
point by a careful finite-size analysis of the order parameter time decay. In particular, 
considering systems' sizes up to L = 4096, we obtained 7 C = 0.32817(2). Being the 
maximum Lyapunov exponent of the single CML A = In 2, from Eq. (j4]) one can derive 
A(7 C ) = —0.3749. Thus the ST takes place at a definitely negative Lyapunov exponent 
and the synchronized state is truly absorbing [7]. 

We evaluate the exponent 9 at criticality by averaging the instantaneous order 
parameter over about 200 independent initial conditions in systems of size L = 1024, 
obtaining a convincing straight line in a doubly logarithmic plot (Fig. Ek)- The best 
fit provides 9 = 0.449(4), which is perfectly compatible with the best DP numerical 
estimates in 2 + 1 dimensions, that is 9 DP = 0.451(6) [5]. The quality of such a scaling 
law is tested by multiplying the order parameter p(t) by t e , as shown in in Fig. [2b we 
obtain an almost two decades long plateau. 

We next compute for several 7 below the critical point at L — 1024. By 
averaging over w 100 different initial conditions we estimate (3 = 0.584(9) (Fig. [2b), to 
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Figure 3. The finite size scaling collapse for ST in CMLs according to Eq. (fTU|) 
is employed to estimate the exponent z. (a) Collapse for CMLs of Bernoulli maps 
obtained using 9 — 0.449 and z = 1.77. Inset: non rescaled data for L = 
16,32,64,128,256 (from left to right), (b) Collapse for skewed tent maps (see Sec. 
I2.2|) obtained using 9 = 1.81, z = 1.55. Inset as in (a). Data have been averaged from 
10 2 up to 10 4 realizations depending on the system size 



be compared with the DP estimate flop = 0.584(4). The larger error is mainly due to 
the uncertainty on the location of the critical point. 

Finally, we determine the dynamical exponent z through the finite-size scaling 
relation ( ITOl . Our best data collapse (shown in Fig. [3k) suggests z = 1.77(3). All 
together, the three critical exponents 9, f3 and z completely identify the DP universality 
class. 

In Tabled] the results are summarized together with the exponent for DP in d = 2 as 
obtained from the best known numerical estimates reported in Ref. [5]. The agreement 
is very good and we can safely affirm that the ST of map with discontinuities in d = 2 
is in the DP universality class, confirming the one dimensional findings. 





9 


P 


z 


Bernoulli 
DP [5] 


0.449(4) 
0.451(6) 


0.584(9) 
0.584(4) 


1.77(3) 
1.76(3) 



Table 1. Critical exponents for coupled Bernoulli maps in 2 + 1 dimensions together 
with the best estimations of the critical indexes for DP in d = 2. For DP 9 and z have 
been derived by employing v\\ = 1.295(8) and v\_ = 0.733(6) [5]. 



2.2. Continuous Maps and Multiplicative Noise universality class 

We now study the synchronization transition in CMLs with continuous local maps. In 
particular, we consider the system ([!]) with local dynamics given by the skewed tent 
map on the unit interval, namely 

ax ifx<l/a 
" S a{x-l)/{l-a) tfx>l/a 1 ' 
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Figure 4. Same as FigQ]for coupled skewed tent maps. The transversal coupling is 
7 = 0.132, slightly above the critical point. The color coding and the scale are the 
same as Fig|T] 

where we set a = 2.2. The skewed tent map is the simplest generic continuous map; 
the skewness ensures the fluctuation of the multipliers (first derivatives of the map) in 
tangent space, which is the generic behavior 6] . Similarly to Fig.[U Fig. H] illustrates the 
spatio-temporal evolution of the synchronization error field wt(x,y) for 7 slightly larger 
than the critical coupling value. Already at first glance, comparing the two figures one 
can argue that the two STs should belong to different universality classes. It is also 
apparent the different nature of the absorbing state (in black). 

To make quantitative the above statement, as for the previous case, we first 
identified the critical coupling 7 C at which ST takes place. 

From the time decay of the order parameter p(t) we estimated 7 C = 0.131760(5). 
This value appears to be compatible with the requirement that the transverse Lyapunov 
exponent At vanishes at the MN synchronization transition. Indeed, from an 
independent numerical simulation of a single replica, we found that the largest Lyapunov 
exponent is A = 0.30588(6), which through Eq. (0} yields A T = 0.00001(6). Thus the 
ST now takes place at a zero transverse Lyapunov exponent. It should be mentioned 
that to minimize finite size effects, which are rather severe in this case (as discussed in 
the following) we have considered lattice up to size L=8192. 

Once the critical point is known, we can estimate the critical exponents. In Fig. Ek,- 
b we report the results on the time decay of the order parameter at criticality. The 
critical exponent 9 is estimated by multiplying p(t) by t 9 and varying 6 as to maximize 
the size of the flat plateau. However, it is worth stressing that, at variance with the case 
of discontinuous maps, here finite-size effects are more severe and numerical artifacts 
[jj may be present. In particular, the asymptotic power law decay sets in at later 
times with respect to discontinuous maps. Moreover, the critical exponent is larger 
than the corresponding DP value, so that lack of statistics tends to plague late time 
data. Therefore, it is necessary to explore large lattice sizes to obtain reasonable scaling 

|| Contrary to what reported in Ref. )23| . we verified that the occurrence of spurious saturation effects 
of the order parameter close to the critical coupling are induced by the employed numerical precision 
and not by long time correlations. 
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Figure 5. Critical behavior of ST for CML's of skewed tent maps (jlip with a = 2.2. 
(a) Power law decay of the order parameter at j c , the dashed line marks the t~ e power 
law with 9 = 1.81. Data have been obtained by averaging over three realizations 
with system size L = 8192 and by a further logarithmic-window average in time, (b) 
The order parameter is multiplied by t to obtain an asymptotically flat curve, (c) 
Subcritical behavior computed from simulations at system size L = 4096. The dashed 
red line marks a slope of 2.19. 



and rule out finite-size effects. We performed numerical simulations in systems of size 
L = 8192, averaging over three independent realizations, obtaining slightly more than 
a decade of convincing scaling behavior. We estimate 9 = 1.81(5). 

As it is shown in Fig. [5b, the behavior of the saturated order parameter in the 
subcritical regime (7 < j c ), allows us to measure the second critical exponent. Our best 
estimates, in systems of linear size L = 4096 give us (3 = 2.19(9). 

For phase transitions in the MN universality class, the dynamical exponent z 
has been conjectured (and confirmed by numerical simulations in 1+1 dimensions) to 
coincide with the one associated to the KPZ equation [35]. Indeed, as it will become 
clearer in Sect. [3], the MN synchronization problem can be mapped in the depinning 
transition of a bounded KPZ surface. From this mapping one deduces that at the critical 
point, where the interface asymptotically depins, MN systems should exhibit the same 
space and time correlations as free KPZ ones. Thus they are also characterized by 
the same z exponent. By rescaling finite size data averaged over many realizations 
with 9 = 1.81 by means of relation ffTUj) . we obtained z = 1.55(8) (see Fig. [3b), 
which is compatible with the known best estimates of KPZ dynamical exponent in 
2+1 dimensions, z KPZ = 1.607(3) as reported in Ref. [36J. 

3. Stochastic Models and Scaling Arguments 

As far as we know, the only measurement of MN critical exponents in two spatial 
dimensions reported in literature [27] was obtained by numerically investigating the 
associated Langevin equation. In order to obtain independent and accurate estimations 
of the critical indexes, we studied two stochastic models which are known, in 1+1 
dimensions, to belong to the MN universality class. 
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3.1. Single Step model plus a hard wall (SSW) 

The MN universality class is closely related to the depinning of a KPZ surface from a 
hard substrate. Indeed the MN Langevin equation [9] 

d t( f> = b(p- c<p 2 + DV 2 <p + (pri (12) 

can be formally mapped, via the Cole-Hopf transformation ip(r,t) = — ln0(r, t), onto 
the KPZ equation with negative nonlinear term and bounded from below [T5] 

dtif) = v - ce~^ + DV 2 ^ - D ( V^) 2 + i] . (13) 

Here, if)(r,t) > and </>(r, t) G [0, 1] are the coarse-grained height and difference field, 
respectively, while T](r,t) is a spatio-temporal Gaussian white noise, and D > 0. Note 
that the </>(r, t) = absorbing state of the Langevin dynamics corresponds to an infinite 
height (i.e. a completely depinned) interface in the bounded KPZ representation. 

The above link suggests us to consider a simple solid-on-solid stochastic deposition 
model, belonging to the KPZ universality class, such as the well known single step model 
(SSM), which in one spatial dimension can be exactly mapped onto the KPZ equation 
(see, e.g., Ref. [37], see also Ref. [38] for a study in two spatial dimensions). Equipped 
with a hard substrate, the so-called SSM-plus-wall (SSW) provides an example of MN 
phase transition. In 1+1 dimension, the critical point is analytically known and very 
accurate numerical estimates of MN critical exponents have been obtained [T4"ll2"6] . The 
SSW time-depinning exponent 9 has also been computed in one spatial dimension via a 
mean field-like approximation in Ref. [39J. 

Here, we numerically investigate the following two dimensional generalization of the 
SSW model. A fluctuating interface with an integer height field h t (x,y) is defined on a 
square lattice (x, y = 1,2, ...,L) with periodic boundary conditions. The dynamics 
of the interface is subjected to the following restrictions: the heights on nearest 
neighbours sites cannot be equal and must differ by one, i.e. \h t (x, y) — h t (x ± l,y)\ = 
\h t (x,y) — h t (x,y ± 1)| = 1, moreover a hard lower wall at height H t moving upward 
with velocity V w is imposed requiring that h t (x,y) > H t , with H t = V w t. The dynamics 
is asynchronous: at each sub-time step dt = 1 /L 2 a site (x, y) is chosen at random and 
its height is increased by two, h t+ dt(x,y) = h t (x,y) + 2, if h t (x,y) is a local minima. 
Finally, the wall moves up by one unit every L 2 /V w time steps pushing up the interface: 
every site whose height is less than H t is raised by two units. The density p of interface 
sites attached to the interface (i.e. h t (x,y) = H t ) immediately after a wall movement is 
the proper order parameter, while the wall velocity V w is the control parameter. 

Unfortunately, the critical velocity V WyC of the two dimensional SSW is not known 
analytically. Careful finite-size analysis up to L = 4096 has been used to locate the 
critical depinning velocity, which is estimated to be V w ^ c = 0.34135(10). Numerical 
simulations start from a completely pinned interface: ho(x,y) = for x + y odd and 
h (x,y) = 1 otherwise. Our results for the time decay critical exponent 9, shown in 
FigEh-b, yield 9 = 1.80(5). By slightly increasing the wall velocity above the critical 
value we can estimate the magnetization exponent, corresponding to the behavior of the 
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Figure 6. Critical behavior of the depinning transition in the d = 2 SSW model 
with L — 4096. (a) Power law decay of the order parameter at criticality, the 
dashed line marks the power law Data are obtained by averaging three different 

realizations and by a further logarithmic-window average in time, (b) Order parameter 
multiplied by t l s to demonstrate the quality of scaling, (c) The subcritical behavior 
is characterized by a power law with exponent (5 = 2.36, as indicated by the dashed 
line. 



ST for subcritical couplings, obtaining (3 = 2.36(9), as shown in Fig [6)3. Summarizing, 
9 is in fairly good agreement with the estimate for the CML with skewed tent maps, 
and (3, while still compatible if error bars are considered, appears to be slightly larger. 
Finally, the dynamical exponent z is estimated via the finite-size scaling relation ffTUl) 
with 9 = 1.8. A satisfactory data collapse can be obtained with z = 1.63(5) (not shown), 
which is compatible with the value obtained for the tent map. 

It is worth concluding by remarking that in Ref. [26J it has been shown that the 
probability distribution of the first depinning time (i.e. the first time at which an 
initially flat and pinned surface depins) does not follow a standard finite scaling relations 
in the MN case. The numerical simulations of the SSW model, however, show that 
other quantities of interest, as the density of pinned sites, follow the typical scaling of 
absorbing phase transitions. This peculiarity can be probably ascribed to the "weakly 
absorbing" nature of the MN absorbing state [33J. 

3.2. Random Multiplier model 

We now consider the Random Multiplier (RM) defined by the dynamics [7] 

aw t (x,y), w.p. ±-aw t (x,y) 

«*»(»,») = < MX ,' y) { A ' f A (15) 

' aw t {x, y), w.p. 1 — aA 

where "w.p." is the shorthand notation for "with probability", w t (x,y) E [0 : 1] and 
w t (x,y) = V^w t (x,y) is the discrete diffusive operator, as in Eq. ([I]), with e = 4/5, on 
a square lattice (x, y — 1, 2, . . . , L) with periodic boundary condition. 
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Before discussing the results of the model, it is worth briefly reviewing its properties 
and its meaning. The model (IT4"1) - (IT5"1) was originally introduced in Ref. [7j to describe 
the synchronization error evolution in proximity of the synchronization transition of both 
continuous and discontinuous maps. Essentially the value of parameter A discriminates 
between discontinuous or continuous character of the local dynamics, switching the 
behavior from the DP to the MN universality class. To better understand the origin of 
the model, first consider Eq. ([1]) equipped with the Bernoulli map ([3]). One can formally 
compute the linearized evolution equation for the synchronization error, which reads 

w t +ii.x,y) = w t (x,y)(l - 2^)d u F(u t (x,y)) . (16) 

Such an equation holds locally for any finite synchronization error wt(x,y) such that 
Ut(x,y) and Vt(x,y) fall on the same branch of the Bernoulli map. In this case one 
simply has w t+ i(x, y) — Wt(x, y) a (1 — 27). However, whenever u t (x, y) and v t (x, y) fall 
on the two different branches of the Bernoulli map, the synchronization error is typically 
expanded to order 1 values, a fact overlooked by the linearization. This latter situation 
can occur with a probability proportional to Wt(x,y) [?]. Setting A = 0, so that only 
Eq. (fT4"j) is relevant, reproduces the Bernoulli map dynamics with a = a(l — 27). 

On the other hand, for the skewed tent map fTTTl) — as well as any continuous map 
- the full dynamics is well captured by linearized dynamics. For the local map ( fTTT) the 
local multiplier assumes one of the two values d u F(a,u) = a or d u F(a,u) = a/(l — a) 
according to the chaotic dynamics. By approximating the chaotic signal with randomly 
chosen multipliers, the RM model with finite A mimics exactly this latter situations, cfr. 
Eq. ( fl5l) ; while Eq. ( fl4l) simply provides a nonlinear saturation effect. Interestingly, in 
Ref. [7J it has been shown that there exist a threshold A > below which the transition 
still belongs to the DP class. This corresponds to the case of almost-discontinuous 
piecewise linear maps on the unit interval, characterized by a very steep branch (with 
slope 1/A and width of order A); for a complete discussion see Ref. |7J and [T2] . 

In the following, being interested to the MN transition, we investigate the RM 
model for A = 0.2, which is a sufficiently large parameter value to drive the system 
into its "linear" regime. The parameter a is the control parameter of the transition; 
indeed, from Eq. ( Tl6l) it follows that a is essentially equivalent to (1 — 2'j)d u F(a, u) so 
that to increase the coupling 7 in the deterministic model amounts to decrease a in the 
stochastic one. The synchronized regime is thus obtained for a smaller than the critical 
value ot c . 

The critical coupling a c implies a zero transverse Lyapunov exponent; it can be 
estimated evaluating the value of a for which (In wt) grows at most logarithmically in 
the linear regime [TJ. This estimation is in agreement with the usual analysis performed 
to evaluate the critical point from the scaling of the order parameter in time (Fig. [7_h). 
In particular, by considering lattices of linear sizes L = 2 x 10 3 and L = 4 x 10 3 the 
best scaling of p{t) is observed at a c = 0.529825(5), in agreement with the Lyapunov 
estimate which gives a c = 0.52981(1). Moreover, the critical exponent is found to be 
6 = 1.76(5), as shown in Fig. [7k-b. For a slightly above a c , from the saturated value of 
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Figure 7. Critical behavior of the RM model in 2+1 dimensions with size L = 2000. 
(a) Power law decay of the order parameter at criticality. Data have been obtained 
by averaging over 100 realizations. The dashed line marks the power law i" 1 - 76 . (b) 
Order parameter multiplied by t 176 ; notice the asymptotic plateau at large times, (c) 
Subcritical behavior computed from up to 400 independent realizations, from which (3 
is estimated to be 2.18 (dashed line). 



the asymptotic density we also find (3 = 2.18(8) (see Fig. [7b). Finally, finite size scaling 
via Eq. ([TO]) gives z = 1.7(1) (not shown). 

3.3. Scaling arguments and comparison among the critical exponents 

This subsection is devoted to a discussion of the results obtained so far for the MN 
universality class. Our three independent estimates of the critical exponents for the MN 
class in 2+1 dimensions are reported in Table [2] For 9 and (3, ST in CMLs, the SSW and 
RM model are in remarkable agreement, the exponents are essentially coincident within 
error bars. In particular, 9 also coincides with the early estimates of Ref . [27] , while a 
remarkable difference (by a factor two) is observed for (3. This is most likely due to the 
preliminary nature of the d = 2 Langevin dynamics simulations [ID] . It should be noted 
that typically a fairly small lack of accuracy in the estimation of the critical point can 
result in a quite large error in evaluating the exponent (3. Also the dynamical exponent 
z measured from finite size scaling essentially agrees with the expected KPZ value, see 
Table El 
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P 


z 




CML (skewed tent maps) 


1.81(5) 


2.19(9) 


1.55(8) 


1.21(8) 


RM 


1.76(5) 


2.18(8) 


1.63(5) 


1.25(8) 


SSW 


1.80(4) 


2.36(9) 


1.7(1) 


1.31(4) 


KPZ [36] and scaling relation [35J 






1.607(3) 


1.32(1) 



Table 2. Critical exponents for CMLs of skewed tent maps, random multipliers model 
(RM) and single step model with a wall (SSW). Also the values that can be deduced 
from scaling arguments and the KPZ numerical estimates are reported. 



Exploiting the relation with free (without walls) KPZ scaling, it has also been 
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conjectured that 

* = W=T) (17) 

which, from zkpz = 1.607(3) as obtained in Ref. [36], implies u± = 0.824(5) and 
z/|| = zu± = 1.32(1). This latter exponent can be compared with the ratio u\\ = (3/9 as 
obtained for the models here analyzed: as reported in Table [2] there is an agreement 
within error bars with the KPZ estimate, as far as the RM and SSW are concerned. 
On the other hand, Uu is slightly underestimated when CMLs of skewed tent maps 
are considered. Altogether, KPZ estimates and the derived (via Eq. (ITTj) temporal 
correlation exponent, well agree with our results in 2 + 1 dimensions for both the 
dynamical exponent z and the ratio (3/9. 



4. Synchronization in the presence of mismatch 

As discussed in the introduction, in typical experimental settings it is impossible to 
produce two exact replicas of the same system. Systematic errors, slight differences 
in the preparation of the system, or inhomogeneous external influence should always 
be taken into account. Such small differences can be mimicked at the level of the 
CMLs model as a quenched random mismatch between the dynamical parameters of 
two otherwise identical systems (this idea was introduced in Ref. [6] in the context of 
low dimensional maps). In practice, we consider the following model, written in generic 
spatial dimension d: 

4%(r) = (l- 1 )F(a 1 (r),u ( t 1 \r))+ 1 F(a 2 (r),u ( t 2 \r)) 

ul%(r) = {l- 1 )F{a 2 {r),u?\r))+ 1 F{a x {r),u?\r)) (18) 

where r is a d dimensional vector and 

fiW( P ) = (i_ e ) u W( r ) + i_ « ( V), (19) 

r'&NN{r} 

where the sum runs over the nearest neighbours of r, denoted as NN{r}. The local 
map F(a,u) depends both on the dynamical variable u = uf\r) and on the quenched 
parameter a = a,(r) (where i — 1,2 labels the replica). Periodic boundary conditions 
are considered as usual in a <i-cube of linear dimension L. We consider two cases, the 
skewed tent map defined in Eq. (ITT]) and the Bernoulli one (Eq. (j3])). Without loss of 
generality, we can write the map parameter a as 

di(r) = a + Ui(r), (20) 

with Ui(r) being a quenched random variable uniformly distributed in [—h, h]. 

Let us now consider the linearized dynamics for the synchronization error w t (r) = 
u\ (r) — uf\r), which reads as (see also [6]) 

w t+1 (r) = w t (r)(l - 2 1 )d u F{a l {r),uf\r)) + (1 - 2 1 )5a(r)d a F(a 1 (r),u ( t 1 \r)) , (21) 

where 5a(r) = ai(r) — a,2(r) is the parameter mismatch. Obviously, ||5a(r)|| ~ h 
(regardless of the chosen norm). The first term on the r.h.s. of Eq. ( 1211) - - i.e. 
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the fluctuating "field derivative" (1 — 2^)d u F{a 1 u) — leads as usual to a stochastic 
term proportional to the amplitude of the synchronization error itself. While the 
second term — i.e. the "parameter derivative" (1 — 2j)d a F(a, u) — , also typically 
fluctuates according to local dynamics, but it depends only on the parameter mismatch 
amplitude and thus acts as an effective "external field" which locally prevent the 
complete synchronization of two replicas. At a field theoretical level, this suggests 
to describe the parameter mismatch as an external driving field with amplitude h. For 
instance, considering the MN class, one can add an additive Gaussian^ white noise ( to 
the Langevin equation (TI2|) 

d t( f) = b(p- c0 2 + £>V 2 + <fyq + h( . (22) 

Implicit in the above formulation is the assumption that the noise terms n and ( are 
completely decorrelated. Of course, the statistical independence between the second 
and the first term in the r.h.s. of Eq. (l2Ti) has to be tested a posteriori. Unfortunately, 
we are not able to solve analytically Eq. (1221) . While efficient numerical methods [4T] are 
known to directly simulate multiplicative noise Langevin equations, we choose instead 
to focus on microscopic models. We consider the RM model (I23H24I) . The external field 
then can be described as an extra additive white noise x{ r ) uniformly distributed in 
[0,h], so that 

( \ J l + Xt(r), w.p. aw t (r) 

wt+xir) = { „ _ if w t [r) > A 23 

aw t {r) + xt{r), w.p. 1 - aw t (r) 

, w t (r) / A + Xt(f) , w.p. aA .„ „ , . . /nd . 

1 awtir) + xt(r), w.p. 1 - aA 

Indeed, the additive noise term x{ r ) creates local activity with an amplitude 
proportional to h, thus preventing complete synchronization between the two replicas. 
In the following we compare direct numerical simulations of the modified RM model 
( 12311241 with the ones of the coupled CMLs (1TB]) with skewed tent maps with a quenched 
mismatch. We expect the asymptotic order parameter p^ to saturate to an /i-dependent 
value, that at criticality should scale as 

Poo ~ h K (25) 

where k is a new critical index (associated to the external field in the field theoretical 
treatment). 

Choosing (aj(r)) = a = 2.2, the critical point estimated in Sec. [2] [j c = 0.13176) 
remains unchanged for h > 0. Equally, we set A = 0.2 for the RM model so that 
a c = 0.52981. 

Starting from random initial conditions, we measure as a function of h in both 
one and two spatial dimensions. The results in two spatial dimensions, reported in 
Fig. [HJi-b, show a remarkable agreement between CMLs and the modified RM model, 

Note that {(} can be always set to zero by an appropriate transformation of the field, <j> ~> 4> + Q- 



Spatio-temporal chaotic synchronization as an absorbing phase transition. 



17 



(a) (b) (c) 
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Figure 8. Critical dependence on the mismatch in CMLs in 2+1 dimensions, (a) 
Mismatch amplitude dependence of the saturated order parameter at criticality for 
skewed maps CMLs (black circles) and RM with an additive noise (red squares). The 
dashed red line marks the h 1 / 2 power law. (b) Same data as (a) rescaled by a factor 
hr x l 2 to highlight critical behavior. System size is L — 2048. (c) For Bernoulli maps 
we check the DP expectation kop = (3/a by plotting p(t)h~ KDP (see text): notice the 
collapse on a single plateau at large times. From top to bottom h = 10~ 7 , 2.5 x 10 -7 , 
6.25 x 10~ 7 , 1CP 6 , 2 x 10~ 6 . System sizes is L = 1024 and data have been averaged 
over 50 initial conditions. Inset: non rescaled data compared with the h = behavior 
(dashed black line). 

thus confirming that the correlations between the additive and the multiplicative terms 
in Eq. ( 12T1) are irrelevant. In particular, we found k = 0.50(5) for the RM and 
k = 0.46(7) for the CMLs. In one spatial dimension (not shown) we have k = 0.21(2) 
(RM) and k = 0.23(2) (CMLs). This is in rough agreement with the d = 1 early 
estimates obtained by direct simulation of the Langevin equation (T22l in Ref. [35J. In 
Ref. [28] a zero dimensional version of Eq. (T22l has been exactly solved via a Fokker- 
Planck approach to yield a logarithmic increase with the field h. While we have not 
been able to put forward any analytical argument in d > 0, these results lead us to 
conjecture k = d/4 for MN in an external field. 

We finally consider the case of Bernoulli map with (aj(r)) = a = 2. Here, we can 
exploit the fact that the scaling theory for Directed Percolation in an external field is well 
known (see [5] and references therein), and thus test if the dependence of on h follows 
the DP prediction. In particular, at criticality, the order parameter should saturates to 
an asymptotic value Poo{h) ~ h K with k = (3 /a, where f3 links to the behavior of the 
density close to the critical point and a controls the mean size of inactive clusters close 
to the critical point, and is related to the other exponents by the formula [5] 

a = u\\ (d/z + 1 — 6) = i/i I + di>± — j3 . 

By means of the best available DP numerical estimates [5], we have k^p = 0.108247(4) 
in d = 1 and Kpp = 0.268(4) in d = 2. As shown in Fig. [Hfc direct numerical simulations 
of 2 + 1 dimensional CMLs with mismatched Bernoulli maps are in very good agreement 
with this prediction. This is also true in 1 + 1 dimensions (not shown). 
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5. Conclusions 

The synchronization transition between spatially extended systems exhibiting spatio- 
temporal chaos is a prototypical example of a fluctuation driven phase transition induced 
by microscopic chaos. The present numerical study in 2+1 dimensions has clearly 
confirmed that, analogously to the 1+1 dimensional case, the synchronization transitions 
can been described in the framework of continuous out-of-equilibrium critical phenomena 
towards an absorbing phase. In particular, depending on the perturbation propagation 
properties of the spatio-temporal dynamics, the synchronization transition belongs to 
two possible universality classes, namely Directed Percolation and Multiplicative Noise. 
The above results confirm that the ST belongs to the DP universality class whenever the 
transverse Lyapunov exponent is negative at the critical point. Differently, MN behavior 
sets in when a zero transverse Lyapunov exponent characterize the critical point. As for 
the latter universality class, we have produced the best available numerical estimates 
of the critical exponents in two spatial dimensions. Furthermore, by analyzing different 
deterministic and stochastic models, we have been able to confirm universality in the 
MN class and its link with deposition processes. 

We have also addressed the effect of a small difference in the dynamics of the 
two replicas, an experimentally relevant question given the practical impossibility to 
produce two exactly identical systems in any experimental setup. By modeling this 
difference as a quenched parametric mismatch in the local dynamics, we have shown 
that it amounts to the action of an external field within the Langevin description. 
Numerical simulations in one and two spatial dimensions for discontinuous local maps 
confirm DP scaling theory, which predicts a power law dependence of the saturated 
density from the external field with an exponent k depending from the usual zero field 
critical ones. Finally, we obtained an estimate for this exponent also for MN in en 
external field (smooth maps). 

Before concluding, a comment on the behavior in higher dimensions is in order. In 
d > 2 naive power counting in the MN Langevin Eq. ( fl2l) predicts the coexistence of 
two different fixed points in the Renormalization Group flow, a mean field one acting at 
small but finite noise amplitude and a strong coupling one at large noise amplitudes [28J. 
Estimates obtained from numerical simulations of Eq. ( fl2l) for MN in 3+1 dimensions 
can be found in Ref. [12]. 

We conclude expressing our hope that these results will stimulate experimental 
studies on the synchronization of extended systems exhibiting spatio-temporal chaos. 
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